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PREFACE 


This  report  presents  major  developments  in  a  continuing  research  program  in  underwater  acous¬ 
tics  at  the  Naval  Research  Laboratory.  These  studies,  begun  in  1974,  have  investigated  the  spatial 
and  temporal  properties  of  low-frequency  acoustic  fields  in  the  ocean.  Their  objective  is  to  provide  a 
priori  estimates  of  the  capabilities  of  highly  directional  receiving  arrays  and  the  signal  coherence  as  it 
may  be  degraded  by  environmental  irregularities.  The  NRL  project  initially  concentrated  on  the 
degradation  caused  by  scattering  in  the  ocean  volume.  These  studies  were  followed  by  studies  of 
scattering  at  its  boundaries.  Both  processes  are  irreversible  and  the  irregularities  causing  the  scatter¬ 
ing  are  statistically  described.  The  statistical  approach  is  natural  for  the  ocean  environment,  because 
the  principal  parameters  exhibit  significant  departures  from  their  mean  values.  They  vary  from  place 
to  place  in  a  manner  that  is  statistically  stable  but  difficult  to  measure  in  detail.  This  is  true  for 
volume  parameters  (e.g.,  sound  speed  deviations)  [1]  and  for  boundary  parameters  (e.g.,  bathymetric 
roughness)  [2],  The  statistical  approach  is  therefore  the  appropriate  method  for  predicting  the  perfor¬ 
mance  of  acoustic  systems  in  the  ocean.  This  is  especially  true  for  receiver  array  design  tasks  such 
as: 

•  Predicting  an  environmentally  constrained  upper  limit  on  array  size  to  achieve  specific 
performance  levels  (gain,  resolution,  and  sidelobe  suppression). 

•  Developing  a  deployment  or  operational  strategy  for  mobile  or  fixed  arrays.  For  example, 
how  many  systems  should  be  used,  where  should  they  be  placed,  and  how  are  they  best 
oriented,  are  all  questions  that  this  statistical  model  can  help  to  answer. 

•  Computing  an  envelope  of  array  system  performance  for  mobile  tactical  systems  where  bot¬ 
tom  contact  at  different  attack  directions  produces  an  output  controlled  by  the  magnitude  and 
orientation  of  bottom  roughness  spectra. 

•  Determining  the  probability  of  detecting  target  echoes  or  passive  radiation  with  specific  sys¬ 
tem  parameters  and  operating  conditions  and  as  a  function  of  the  bottom  spectral  character. 

The  NRL  project  has  addressed  the  causes  of  degradation  of  spatial  coherence  of  acoustic  signals 
as  the  key  new  system  issue.  This  is  because  substantial  investment  in  technology  depends  on 
whether  acoustic  field  properties,  driven  by  tactics,  systems,  and  operating  area,  will  allow  the 
coherent  gain  that  antennas  must  achieve  to  be  effective.  Ocean  acoustic  volume  effects  have  been 
dealt  with  in  earlier  reports  [3-5].  This  research  should  be  extended,  because  scattering  mechanisms 
and  higher  frequency  regimes  that  were  not  previously  studied  are  now  of  significant  interest. 
Although  ocean  volume  scattering  is  always  present  in  underwater  acoustic  propagation,  it  is  generally 
less  than  bottom  scattering  at  all  frequencies.  The  detailed  bottom  irregularities  affecting  an  acoustic 
signal  at  any  given  instant  usually  cannot  be  known  in  sufficient  detail  to  correct  for  the  distortion 
they  cause,  but  it  is  important  to  reliably  estimate  the  expected  reduction  in  receiver  array  perfor¬ 
mance.  This  report  presents  a  model  whose  current  version  provides  estimates  of  beam  spreading  of 
array  outputs  caused  by  a  single  encounter  with  a  hard,  rough  ocean  bottom. 

This  report  describes  the  theoretical,  numerical,  and  computational  basis  of  a  computer  model 
called  the  ocean-refraction  bathymetric-scattering  (ORBS)  model.  This  model  estimates  the  effect  of 
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scattering  by  a  rough,  stochastically  characterized,  ocean  bottom  on  the  three-dimensional  spatial 
coherence  of  an  acoustic  signal.  The  model  allows  for  a  user-defined  sound-speed  profile  in  the 
ocean.  This  model  allows  refraction  in  the  water,  three-dimensional  bathymetric  scattering,  and  spec¬ 
ular  reflection  to  be  combined  realistically  in  the  calculation.  This  combination  gives  ORBS  a  unique 
computational  predictive  capability  that  was  previously  unavailable.  The  model  displays  the  computed 
acoustic  field  at  a  specified  receiver  location  in  terms  of  the  simulated  output  of  an  arbitrarily  oriented 
line  array  centered  there.  (With  minor  modifications,  the  program  could  model  a  more  general  array 
type,  such  as  a  line  array  with  missing  elements  or  a  billboard  array.)  The  computed  outputs  include 
the  main  peak  and  sidelobe  levels  of  the  simulated  pattern,  the  half-power  width,  and  the  computed 
source  bearing  relative  to  the  array. 

In  its  present  form,  the  model  considers  cases  where  one  bottom  interaction  is  dominant.  With 
minor  modifications,  the  ORBS  model  can  be  adapted  to  some  straightforward  variations  of  this  prob¬ 
lem,  such  as  single-bounce  propagation  under  rough  ocean  ice. 
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THE  OCEAN-REFRACTION  BATHYMETRIC-SCATTERING  (ORBS)  MODEL 


1.  INTRODUCTION 

When  a  sound  wave  impinges  upon  a  hard  rough  surface,  a  significant  amount  of  the  incident 
energy  is  incoherently  scattered.  This  is  an  important  phenomenon  that  degrades  the  performance  of 
acoustic  receiving  systems  operating  on  a  signal  field  dominated  by  reflection  and  scattering  from 
rough  bathymetry.  An  array  located  over  a  rough  bottom  may  be  insonified  by  direct  refracted 
arrivals,  by  reflections  from  areas  on  the  bottom,  and  by  energy  scattered  from  the  bottom.  Gen¬ 
erally,  depending  on  the  acoustic  frequency,  the  source/receiver  geometry,  and  the  characteristics  of 
the  bottom,  any  one  of  these  components  of  the  received  sound  may  be  dominant.  Consequently,  if 
the  degradation  in  array  performance  caused  by  bottom  scattering,  which  is  usually  a  very  strong 
diffusing  mechanism,  is  to  be  evaluated,  the  contribution  of  each  component  must  be  computed  accu¬ 
rately. 

The  ocean-refraction  bathymetric-scattering  (ORBS)  model  computes  array  response  in  the 
acoustic  field  by  propagation  of  the  Fourier-transformed  two-point  mutual  coherence  function  of  the 
pressure  field  (the  transform  can  be  thought  of  as  the  power  distribution  in  location  and  direction) 
along  its  characteristic  trajectories.  To  compute  these  refracted  trajectories  accurately,  ORBS  allows 
for  a  realistic  depth-stratified  sound-speed  profile  as  input  data;  a  continuous,  stratified  curvilinear 
interpolating  function  then  integrates  the  path  elements.  The  number  of  paths  required  is  dictated  by 
the  number  of  receiver  points  on  the  computed  array,  the  number  of  beam  directions  considered,  and 
the  need  to  sample  the  acoustic  field  finely  enough  to  integrate  the  beam  power  accurately.  In  a  typi¬ 
cal  problem,  tens  or  hundreds  of  thousands  of  characteristics  may  be  required.  For  each  scattered 
path,  a  cross  section  (probability  for  scattering)  must  be  computed;  ORBS  uses  a  modified  form  of 
Kirchhoff  theory  [6]  for  this.  To  make  the  large  number  of  cross-section  calculations  feasible,  the 
model  uses  an  approximation  to  the  Kirchhoff  integral  that  replaces  the  numerical  integration  by  a  sin¬ 
gle  function  evaluation  (of  the  bottom  roughness  spectrum).  This  reduces  the  computation  time  by 
almost  two  orders  of  magnitude  with  little  loss  of  accuracy  [7]. 

Because  it  must  deal  with  scattering  out  of  the  plane  of  the  incident  wave,  the  calculation  of  the 
characteristics  is  three-dimensional.  To  ‘start  up’  this  three-dimensional  calculation,  ORBS  ordinarily 
uses  a  parabolic  equation  (PE)  model  [8]  to  propagate  the  sound  from  the  source  to  the  scattering 
region.  In  a  long-range  case,  propagation  up  to  the  region  of  the  receiver  is  typically  modeled  as 
cylindrically  symmetric,  and  a  split-step  PE  calculation  in  the  range-depth  plane  is  used  [9,10],  How¬ 
ever,  three-dimensional  PE  models  exist  [11]  and  could  be  used,  if  appropriate,  for  this  phase  of  the 
problem.  In  a  typical  short-range,  bottom-bounce,  active  or  passive  sonar  case,  the  startup  field  is 
required  at  a  range  of  a  few  kilometers  (or  less)  from  the  source,  and  a  cylindrically  symmetric  PE 
calculation  or  simulated  source  field  is  used. 
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Section  2  is  an  overview  of  the  propagation  of  the  coherence  function  transform  in  ORBS  and  its 
use  in  computing  the  directional  response  of  an  array.  Section  3  describes  the  model  formulation  for 
the  redirection  of  acoustic  energy  at  a  scattering  point  on  the  bottom.  Section  4  presents  sample  com¬ 
puted  results  for  varying  geometries  and  varying  parameters  of  the  two-dimensional  spectrum  describ¬ 
ing  the  rough  bottom.  Section  5  is  a  summary  of  the  report.  The  important  computational  algorithms 
used  are  described  in  detail  in  appendixes,  and  a  guide  to  program  operation  is  included. 

2.  THE  MUTUAL  COHERENCE  FUNCTION  AND  ARRAY  RESPONSE 

Figure  1  shows  the  two  phases  of  the  ORBS  propagation  calculation.  The  first  phase  comprises 
the  region  from  the  source  to  the  range  Rp  in  the  figure.  In  a  short-range  bottom-bounce  case  this 
region  may  span  only  a  few  kilometers  and  serve  simply  to  start  up  the  bottom-interaction  phase  of 
the  calculation.  In  a  long-range,  deep-basin  case  where  a  sound  channel  is  present,  it  may  span  tens 
or  hundreds  of  miles.  The  main  phase  of  the  calculation  is  in  the  region  where  bottom  interaction— in 
particular,  a  single  bottom  bounce— is  expected  to  provide  the  dominant  paths  for  sound  reaching  the 
receiver.  Bottom  reflections  and  scattering  introduce  propagation  directions  that  are  not  cylindrically 
symmetric  about  the  source,  so  a  fully  three-dimensional  calculation  is  required  in  this  region. 
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Fig.  1  —  The  concept  of  the  ORBS  propagation  algorithm  in  the  bottom 
interaction  region.  Rp  is  the  “source  cylinder”  range,  and  Rc  is  the 
range  of  the  center  of  the  array.  A  scattered  path  and  a  corresponding 
set  of  incident  paths  are  shown. 


In  the  first  region,  before  significant  bottom  interaction  is  encountered,  a  two-dimensional  model 
is  usually  adequate.  To  start  up  the  ORBS  algorithm,  the  two-dimensional  Helmoholtz  equation  for 
the  complex  acoustic  pressure  p(r,  z)  is  approximated  in  this  region  by  the  parabolic  equation  (PE) 
[8] 


2ik0pr  +  Px  +  H[n\r,  z)  -  l]  p  =0.  (1) 

The  coordinate  system  consists  of  depth  z  measured  positive  downward  from  the  ocean  surface  and 
range  r  measured  horizontally  from  the  sound  source.  In  Eq.  (1),  n  is  the  index  of  refraction 
(n(r,  z)  =  c0/c(r,  z),  where  c(r,  z)  is  the  sound  speed  at  the  point  r,  z  and  c0  is  a  reference 
sound  speed),  and  k0  =  2tt/ /c{)  where  /  is  the  source  frequency  measured  in  Hz.  The  parabolic 
equation  is  solved  by  the  split-step  algorithm  [9,10]. 
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Thus,  with  use  of  the  parabolic  equation  method,  the  acoustic  field  is  known  at  range  r  —  Rp 
(Fig.  1)  as  a  function  of  depth.  At  this  range,  we  form  the  two-point  product  of  the  acoustic  field, 

r  ( r ,  Zi,  z2)  =  p(r,  Z\)p *  (r,  z2), 

where  the  asterisk  indicates  complex  conjugation.  The  two  points  are  at  the  same  range  r .  The  prod¬ 
uct  is  now  written  in  mean-and-difference  coordinates  in  depth: 

z  =  (z i  +  z2)/ 2,  s  =  z\  -  z2 

and  transformed  in  the  depth-difference  coordinate 


r  (r,  z,  <t>)  =  f  r(r,  Z  +  2  -  h  e  isk^ds  . 

•'-»  2  2 


(2) 


T(r,  Z\,  z 2)  is  the  two-point  mutual  coherence  function  of  the  complex  acoustic  pressure  field. 
Its  Fourier  transform  r(r,  z,  <t>)  can  be  regarded  as  a  power  density  or  distribution  in  <j>,  with  o  being 
the  direction  of  energy  propagation  in  the  vertical.  This  function  is  useful  because  there  is  a  transport 
equation  for  propagating  it,  and  there  is  a  method  of  using  it  to  compute  array  response.  In  the  text 
that  follows,  the  transformed  coherence  function  is  denoted  as  f(r,k)  where  r  denotes  a  point  in 
three-dimensional  space  and  the  unit  vector  k  specifies  a  direction.  The  relation  between  f  in  three 
dimensions  and  I  defined  in  a  plane  is  discussed  rigorously  in  Appendix  A. 

From  the  Helmholtz  or  the  parabolic  equation,  a  transport  equation  can  be  derived  for  propagat¬ 
ing  the  transformed  coherence  function.  This  is  a  differential  equation  whose  solutions  are  the 
characteristic  contours,  along  each  of  which  the  value  of  f  is  constant.  References  5  and  12  provide 
further  discussion  of  the  transformed  coherence  function  and  its  transport  equation. 

In  the  bottom-interaction  region,  ORBS  uses  a  propagation  algorithm  based  on  the  fact  that  the 
characteristics  of  the  transport  equation  for  f  are  trajectories  computationally  identical  to  the  rays  of 
geometric  acoustics  [12].  Within  this  limited  region,  we  assume  that  horizontal  variation  of  the  sound 
speed  in  the  ocean  is  negligible,  i.e.,  the  sound  speed  is  a  function  of  depth  only.  The  angular  power 
distribution  at  a  point  on  the  array  is  given  by  the  value  of  f  at  the  point.  This  is  computed  by  trac¬ 
ing  a  discrete  set  of  characteristics  (rays)  from  points  along  the  array  backward  to  the  “source 
cylinder”  where  the  characteristic  value  of  f  for  each  path  can  be  evaluated.  This  source  cylinder  is 
defined  by  r  —  Rp .  To  compute  the  total  field,  it  is  necessary  to  include  contributions  from  paths 
specularly  reflected  from  the  bottom,  from  paths  scattered  by  it,  and  from  direct  waterborne  paths 
(those  not  encountering  the  bottom).  The  backward-looking  algorithm  is  used  for  computational  effi¬ 
ciency;  it  guarantees  that  the  program  never  needlessly  computes  paths  that  do  not  reach  the  array. 

The  program  thus  computes  f1  at  a  set  of  discrete  points,  which  can  be  regarded  as  hydrophones 
in  an  array.  (We  consider  a  linear  array  here  for  simplicity  and  because  it  represents  the  current 
implementation  of  the  model.)  The  next  step  is  to  compute  the  magnitude  of  the  directional  array 
response,  usually  defined  as 

L/  2 

A(d)  =  |  w(Y)  p(Y)  exp  (-ik0Y  sin  6)  dY ,  (3) 

-L/2 
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where  Y  is  the  coordinate  locating  a  point  on  the  array,  p(Y )  is  the  acoustic  pressure,  k0  is  the  acous¬ 
tic  wave  number,  d  is  the  steering  angle,  L  is  the  array  length,  and  w(Y)  is  the  weighting  function 
used  in  the  beamforming.  Appendix  B  shows  that  the  angular  power  distribution  observed  by  the 
array  can  be  written  as  a  weighted  integral  of  F  over  the  array: 

<  \A(d)\2>  =  \  j  j  a(Y,  6)  f  (Y,  k)  cos  ^  dU  dY ,  (4) 


where  a(Y ,  6)  is  computed  from  w(y).  This  observed  intensity  distribution  is  an  ensemble  aver¬ 
age  arising  from  the  stochastic  description  of  the  rough  bottom.  The  angles  \ p  and  Q  specify  the 
direction  of  k  in  a  coordinate  system  to  be  defined  next. 

The  array  may  be  yawed  and/or  tilted  with  respect  to  the  fixed  coordinate  system  in  the  ocean. 
It  is  convenient  to  define  directions  of  paths  traced  from  a  point  on  the  array  in  a  coordinate  system 
that  is  fixed  with  respect  to  the  array.  Figure  2  shows  the  unit  vector  k',  the  initial  direction  of  a 
path  traced  from  point  P  on  the  array.  (The  unit  vector  k  =  — k'  is  the  local  direction  of  the  wave 
vector  on  this  path.)  In  Fig.  2,  ^  is  measured  from  a  plane  whose  normal  is  the  array  axis;  it  is  a 
beam  angle  measured  from  broadside.  All  unit  vectors  k'  having  a  common  beam  angle  describe  a 
right  circular  cone,  and  Q  is  the  angle  defining  the  direction  within  the  cone.  These  coordinates  are 
chosen  because  a  linear  array  cannot  distinguish  among  different  arrivals  having  a  common  beam 
angle.  ORBS  samples  the  acoustic  field  by  tracing  out  a  two-dimensional  array  of  paths  at  equal 
increments  in  \j/  and  fi.  Each  path  is  traced  to  its  intersection  with  the  bottom  (paths  not  interacting 
with  the  bottom  are  a  special  case  to  be  discussed  below).  Appendix  C  describes  the  tracing  formulas 
used  in  the  model. 


Fig.  2  —  Angular  coordinates  defining  a  path 
direction  with  respect  to  the  array.  All  unit 
vectors  kr  at  point  R  on  the  array,  with  beam 
angle  form  a  cone.  Contributions  from  all 
angles  12  within  the  cone  are  summed. 


Because  ORBS  uses  a  receiver-to-source  propagation  algorithm,  it  expresses  the  power  prop¬ 
agated  by  each  scattered  path  as  a  weighted  summation  over  all  possible  paths  from  the  source  to  the 
scattering  point.  The  weighting  function  in  this  summation  is  the  differential  cross  section  for  scatter¬ 
ing  from  the  rough  bottom  and  is  discussed  in  Section  3.  The  program  also  tests  for  the  existence  of 
a  specular  reflection  at  each  bottom  scattering  point  and  computes  the  appropriate  reflection  coeffi¬ 
cient  whenever  a  reflected  path  is  found. 

The  special  case  of  paths  not  interacting  with  the  bottom  is  simpler  because  each  path  can  be 
traced  directly  back  to  the  source  cylinder  in  a  single  vertical  plane.  (This  follows  from  the  assump¬ 
tion  that  the  sound  speed  depends  only  on  depth  in  the  bottom-bounce  region.)  The  refracted-only 
part  of  the  calculation  accounts  for  a  small  part  of  the  computer  time,  most  of  which  goes  to  comput¬ 
ing  the  large  set  of  incident  paths  for  each  scattered  path. 
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The  essential  features  of  the  ORBS  propagation  algorithm  can  be  summarized  as: 

•  A  wave  model  is  used  to  propagate  the  acoustic  field  from  the  source  to  the  bottom-bounce 
region.  At  this  point,  the  mutual  coherence  function  and  its  Fourier  transform  F  are  com¬ 
puted,  because  a  raylike  transport  equation  for  f  exists  and  because  f  can  be  used  to  compute 
array  response. 

•  Ray  computations  are  used  to  propagate  the  field  through  this  bottom-bounce  region,  so  that  a 
fully  three-dimensional  calculation  of  the  refracted,  coherently  reflected,  and  incoherently 
scattered  paths  can  be  made  there.  Also,  the  array  and  the  mean  bottom  plane  may  both  be 
oriented  arbitrarily.  In  principle,  any  array  configuration  could  be  represented  by  a  straight¬ 
forward  modification  of  the  program. 

•  The  raylike  computation  of  bottom-scattered  paths  permits  the  use  of  Kirchhoff  theory  to  com¬ 
pute  the  scattering  cross  section  as  a  weight  for  each  path  when  the  paths  are  summed  at  the 
array.  Thus  the  resulting  summation  includes  a  stochastic  ensemble  average  of  the  angular 
power  distribution  of  bottom-scattered  energy  received  at  the  array.  This  distribution  function 
is  used  to  compute  a  simulated  array  response. 

3.  COMPUTING  ARRAY  RESPONSE  WITH  BOTTOM-INTERACTING  PATHS 

This  section  shows  how  the  array  response  integral  of  Eq.  (4)  is  formulated  when  the  function  f 
in  the  integrand  is  computed  from  paths  that  are  bottom-reflected  or  scattered  between  the  source  and 
the  array.  These  calculations  require  transformations  between  different  coordinate  systems  and 
transformations  of  the  arguments  of  delta  functions  in  the  integrands,  which  may  introduce  nontrivial 
Jacobian  factors. 

Figure  3  shows  the  geometry  of  the  problem.  We  assume  that  a  single  sound  speed  profile 
characterizes  the  bottom-interaction  region,  so  that  each  path  traced  (at  least  until  it  contacts  a  boun¬ 
dary)  is  confined  to  one  vertical  plane.  This  plane,  labeled  PR,  contains  a  point  R  on  the  array  and  a 
scattering  point  B  in  the  mean  bottom  plane.  (The  angle  a  is  the  yaw  of  the  array.)  The  vertical 
plane  Ps  contains  B  and  the  source  5.  The  angles  6  and  8  specify  the  azimuths  of  the  two  planes, 
and  nft  and  n5  are  the  respective  unit  normals.  Thus  the  figure  defines  the  propagation  planes  for 
two  traces,  Ps  for  an  incident  path  and  PR  for  a  scattered  (or  reflected)  path.  The  two  paths  meet  at 
point  B,  as  shown  in  Fig.  4.  Unit  vectors  q  and  k  are  the  directions  of  the  incident  and  scattered 
paths,  respectively,  at  B . 

It  is  convenient  to  define  the  directions  of  q  and  k  in  a  coordinate  system  that  is  fixed  with 
respect  to  the  bottom  at  B ,  because  the  scattering  cross  section  is  a  function  of  the  incident  and  scat¬ 
tered  angles  in  such  a  system.  Figure  4  shows  the  system  we  use;  its  basis  is  the  set  of  unit  vectors 
[n,  p,  d],  where  n  is  the  normal  to  the  mean  bottom  plane,  p  is  the  downslope  direction  in  the  mean 
bottom  plane,  and  d  =  n  x  p  (making  d  horizontal).  The  incoherent  part  of  f  in  the  scattered 
direction  can  be  expressed  as  an  integral  over  incident  directions,  where  6  is  the  polar  angle  meas¬ 
ured  from  d  and  <f>q  is  an  azimuth  measured  from  p: 

7T 

Tsca,  <rB ,  k)  =  fQ  sin  dq  |  ^  k)  S(q-n,)f(rB,  q)  d<j>q d9q  ,  (5) 

2 

where  rB  is  the  location  of  the  scattering  point,  a  (q  —  k)  is  the  differential  cross  section  per  unit 
area  for  scattering  from  q  into  k  (see  Appendix  D),  and  the  delta  function  requires  that  each  incident 
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Fig.  3  —  The  coordinate  system  used  in  the  propa¬ 
gation  calculation.  The  £-axis  is  positive  downward 
(into  the  page).  The  angles  <5  and  Q  are  azimuths  of 
the  vertical  planes  of  the  antecedent  and  scattered 
paths,  respectively.  The  angle  a  is  the  yaw  of  the 
array  (rotation  about  a  vertical  axis,  from  broad¬ 
side). 


Fig.  4  —  The  antecedent  and  scattered  directions  at 
the  bottom  point  B  are  defined  by  wave  vectors  q 
and  k,  respectively.  The  bottom-fixed  coordinate 
system  [n,  p,  d]  is  shown. 


vector  q  must  lie  in  Ps.  The  integral  of  Eq.  (5)  can  be  evaluated  by  some  manipulations  that  simplify 
the  delta  function  so  that  the  integration  over  <j>q  can  be  performed.  The  resulting  integral  over  0  is 
easily  transformed  to  an  integral  over  <f>,  the  vertical  angle  of  the  incident  direction  in  Ps : 

7 r 

fscat  (rB ,  k)  =  j  2  -g(q  f  k)  f  (rB ,  qm  d<p.  (6) 

n  ■  k 

2 

From  Eq.  (4)  the  array  response  as  a  function  of  steering  angle  6  caused  by  the  incoherent  part  of  the 
acoustic  field  is  written 

<  M(0)|  ^  scat  HI  a(r,  \p,  6 )  f  scat  (rB(y ,t/sO),  k  ( Y,\J/,Q )  cos  \p  d\p  dQ  dy  (7) 

where  Tscat  (rB,  k)  is  given  by  Eq.  (6).  Equation  (7)  explicitly  shows  the  dependence  of  rB,  q,  and  k 
on  Y,  \p,  and  0  (the  location  and  direction  of  the  scattered  path  at  the  array). 

Now  consider  the  contribution  of  the  specularly  reflected  part  of  the  field,  which  can  be  written 
(from  Eq.  (4)) 


</l(0)|2>refl  =  j  jJfl(y,Wf  (rB(y,^Q),  k(y,^,fi))  cos  td+dtldY.  (8) 
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A  A 

When  T(rB ,  k),  the  characteristic  value  on  the  reflected  path,  is  computed  by  evaluating  the  value  on 
the  incident  path,  Eq.  (8)  becomes 

<  I A  (0)  | 2 > ref]  =  j  J  J  a(Y,t,d)  f  k (E,^,Q))  5  (qrefl  -  ft,)  cos  ^d^dSldY  ,  (9) 


where  qrefl  must  be  in  the  specular  direction: 

qrefl  (y,^,fi)  =  k (y,^,Q)  -  2n  [n-  k(Y,t,Q)]. 


Again,  the  notation  indicates  that  the  location  rB  and  the  directions  k  and  q  of  the  reflected  and 
incident  paths  at  B  are  functions  of  the  location  and  direction  of  the  reflected  path  at  the  array.  The 
delta  function  requires  that  the  incident  path  producing  the  reflection  be  in  the  vertical  plane  contain¬ 
ing  the  source;  otherwise  the  integral  must  vanish.  As  above,  some  manipulation  of  the  delta  function 
is  used  to  facilitate  the  integration  over  one  variable  in  the  multiple  integral.  The  result  is  an  integral 
over  fi  and  Y : 


<  M(0)l2> 


r(rB(y,^Q),Q),q(y,^(fi),«) 


^•*V3refl(y,» 


cos  dQ  dY . 


(10) 


Here,  the  derivative  in  the  denominator  can  be  approximated  numerically  by  tracing  a  path  having  a 
slightly  different  value  of  i p  whenever  a  specular  reflection  is  found.  Ordinarily  only  one  solution 
\pi  (U)  exists. 

4.  SAMPLE  OUTPUTS  FROM  THE  ORBS  MODEL 

In  this  section  we  present  the  results  of  some  sample  ORBS  calculations.  The  first  case  uses 
environmental  data  from  the  vicinity  of  the  Gorda  Rise  in  the  northeastern  Pacific  (near  43  °N, 
127°W),  an  area  whose  bathymetry  has  been  extensively  mapped  and  statistically  analyzed  [2,13], 
Roughness  “provinces”  in  the  region  have  been  found  with  spatial  amplitude  spectra  having  a 
power-law  dependence,  in  the  formulation  of  Ref.  2, 

A(k)  ~  k~b  , 


where  k  is  the  spatial  wave  number  and  1  <  b  <  2.  An  amplitude  spectrum  in  Ref.  2  with  a  power 
law  coefficient  b  =  -1.5  corresponds  to  a  roughness  spectral  density  function  in  our  formulation 
(see  Appendix  D) 


S(k)  ~  k~4. 


Accordingly,  we  model  the  ocean  bottom  in  a  typical  subarea  of  the  Gorda  Rise  with  the  roughness 
spectrum  and  correlation  function  (its  two-dimensional  Fourier  transform) 

S(k )  =  47r(aL)2/[l  +  (kL  )2]2 , 
and 


W(r)  =  (r/L)  K{(r/L), 
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where  a  is  the  rms  residual  bottom  height  and  L  is  a  correlation  length.  L  may  be  taken  to 
correspond  to  the  longest  spatial  wave  length  observed  in  the  subarea  (i.e.,  its  size)  and  a  is  the  total 
rms  residual  height  over  all  A: -space.  is  the  modified  Bessel  function  of  the  second  kind. 

Figure  5  shows  the  geometry  used  for  the  sample  calculations.  The  mean  ocean  bottom  is  a 
horizontal  plane  of  depth  2972  m.  A  1  kHz  beamed  source  5  is  at  a  depth  of  100  m.  The  receiving 
array  is  at  C,  also  at  a  depth  of  100  m.  The  source  beam  pattern  is  aimed  15°  down,  with  a  half¬ 
power  width  of  5°.  Figure  6  is  a  gray-scale  plot  showing  the  result  of  a  parabolic-equation  intensity 
calculation  for  this  problem,  with  a  highly  reflective  bottom  (20°  critical  angle).  The  shading  in  each 
range-depth  cell  denotes  the  intensity  relative  to  cylindrical  spreading;  darker  shading  means  higher 
intensity.  Figure  7  is  a  similar  plot.  However,  for  the  PE  calculation,  bottom  roughness  here  is 
introduced  by  pseudorandom  variation  of  the  interface  depth  about  the  mean  bottom  depth,  with  an 
rms  residual  roughness  of  20  m  and  a  correlation  length  of  400  m.  The  rough  bottom  thus  generated 
is  one  realization  of  an  ensemble  of  bottoms  characterized  by  these  statistical  parameters.  (The  ORBS 
calculations  to  follow  are,  of  course,  based  upon  the  statistics  of  ensemble  averages.)  Although  the 
two  figures  have  the  same  general  features,  Fig.  7  clearly  shows  bottom  scattering  and  qualitatively 
suggests  deterioration  of  coherence  in  the  reflected  and  scattered  field. 


1000-Hz  Source  150-m  Array 


0  20  km 


Fig.  5  —  Configuration  for  the  sample  calculations 
of  Figs.  6  to  16.  A  beamed  1  kHz  source  is  at  5, 
and  the  receiving  array,  which  may  be  horizontal  or 
vertical,  is  centered  at  C. 


3200  1 - : — - — - : - - - J 

0  5  10  15  20  25 


Range  (km) 

Fig.  6  —  Gray-scale  intensity  plot  for  the  sample 
problem,  from  a  parabolic  equation  calculation  with  the 
Gorda  Rise  environment  and  a  smooth  reflective  bottom 


3200  - - - ■ - - - — i— 

0  5  10  15  20  25 

Range  (km) 

Fig.  7  —  Gray-scale  intensity  plot,  similar  to  Fig.  6, 
with  the  smooth  bottom  replaced  by  a  deterministic  rough 
bottom 
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The  receiver  in  the  problem  is  a  150-m  array  centered  at  a  depth  of  100  m  and  oriented  either 
vertically  or  horizontally.  In  Figs.  8  and  9  the  array  is  at  a  range  of  20  km,  in  the  beam  of  bottom- 
reflected  energy,  and  horizontal.  In  Fig.  8  ORBS  has  averaged,  over  the  21  phones  in  the  array,  the 
relative  beam  levels  for  the  smooth-bottom  case  of  Fig.  6,  i.e.,  the  received  energy  reaches  the  array 
by  coherent  reflection.  Figure  9  is  an  analogous  plot  for  a  rough-bottom  case  with  the  parameters  of 
Fig.  7,  i.e.,  the  received  energy  reaches  the  array  by  bottom-scattered  paths.  (The  power  levels  on 
these  plots  are  relative  to  the  same  arbitrarily  normalized  source  level.)  The  peak  of  Fig.  9  is  ~8  dB 
lower  and  2  broader  than  the  peak  of  Fig.  8.  Figure  10  shows  the  array-processed  beam  patterns 
corresponding  to  the  two  cases  superimposed  on  a  single  graph.  (Appendix  B  describes  our  algorithm 
for  array  simulation.)  The  half-power  width  increases  from  0.4°  for  the  reflected  case  to  2.5°  for  the 
scattered  case. 
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Fig.  8  Directional  distribution  of  received  energy  with 
horizontal  receiving  array  at  20-km  range  and  smooth 
reflective  bottom 


In  Figs.  11  and  12,  the  array  is  vertical  and  the  mean  phone  levels  are  computed  (as  in  Figs  6 
and  7).  Figure  11  is  the  smooth-bottom  case  at  20  km,*  and  Fig.  12  is  the  rough-bottom  case.  Once 
again,  the  peak  is  lowered  and  broadened  in  the  rough-bottom  case.  Figures  13  and  14  show  the  hor¬ 
izontal  array  at  a  range  of  15  km,  at  the  edge  of  the  zone  between  the  source  and  the  first  strong 
bottom-reflected  returns.  Not  surprisingly,  the  received  level  in  the  reflected  case  (Fig.  13)  is  8  dB 
lower  than  at  20  km.  However,  in  the  scattered  case  (Fig.  14),  the  peak  level  has  dropped  only  3.6 
dB,  showing  that  a  rough  bottom  scatters  a  significant  amount  of  energy  into  the  otherwise  uninsoni- 
fied  region.  In  fact,  at  shorter  ranges  in  this  region,  scattered  energy  is  received  where  the  reflected 
energy  has  virtually  disappeared.  Figure  15  summarizes  this  effect  by  plotting  the  relative  maximum 
beam  levels  at  the  horizontal  array  as  functions  of  range  for  the  reflected  and  scattered  cases.  The 
reflected  energy  disappears  for  ranges  <  14  km,  while  in  the  scattered  case  energy  still  appears  here. 


*The  sharp  cutoff  at  17°  in  Fig.  11  results  from  the  cutoff  in  the  eigenray  set  in  the  smooth-bottom  case. 
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Fig.  9  —  Directional  distribution  of  received  energy  for  the 
sample  problem  with  the  horizontal  receiving  array  at  20  km 
range  and  stochastic  bottom  scattering.  The  source  level  is 
the  same  as  in  Fig.  8. 


Beam  Angle  (deg) 

Fig.  10  —  Array-simulated  beam  response  patterns 

corresponding  to  Figs.  8  and  9.  The  dashed  curve  is  for  the 
smooth  bottom  (a  =  0)  and  the  solid  curve  is  for  the  rough 
bottom  (a  —  20  m). 
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Fig.  11  —  Directional  distribution  of  received  energy  for  the 
sample  problem  with  vertical  receiving  array  at  20  km  range 
and  smooth  reflective  bottom.  The  source  level  is  the  same  as 
in  Fig.  8. 


Fig.  12  —  Directional  distribution  of  received  energy  for  the 
sample  problem  with  vertical  receiving  array  at  20  km  range 
and  stochastic  bottom  scattering.  The  source  level  is  the  same 
as  in  Fig.  8. 
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Fig.  13  —  Directional  distribution  of  received  energy  for 
the  sample  problem  with  vertical  receiving  array  at  15 
km  range  and  smooth  reflective  bottom.  The  source 
level  is  the  same  as  in  Fig.  8. 


1000  Hz  150-m  array  at  100  m 

Array  tilt  0.0°,  Yaw  0.0° 

Bottom  slope  0.0°,  Yaw  0.0° 
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Fig.  14  —  Directional  distribution  of  received  energy  for 
the  sample  problem  with  vertical  receiving  array  at  15 
km  range  and  stochastic  bottom  scattering.  The  source 
level  is  the  same  as  in  Fig.  8. 
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Fig.  15  —  Level  of  peak  arrival  at  horizontal  array 
as  a  function  of  range  for  a  smooth  bottom  (a  =  0) 
and  a  stochastic  scattering  bottom  ( o  =  20  m)  in  the 
shadow-zone  region  for  the  Gorda  Rise  case 


Figure  16  shows  the  geometry  for  a  second  set  of  calculations.  A  100  Hz  source  is  at  a  depth 
of  100  m,  216  km  from  a  360-m  horizontal  array  at  a  depth  of  2700  m.  The  bottom  slopes  down¬ 
ward,  away  from  the  array,  at  an  angle  of  15°.  The  bottom  is  characterized  by  the  same  isotropic 
correlation  function  used  in  the  short-range  problem,  but  both  a  and  L  are  allowed  to  vary.  Figure 
17  shows  a  family  of  plots  of  array  signal  gain  (ASG)  vs  rms  roughness,  with  the  bottom  correlation 
length  as  a  parameter.  Figure  18  is  an  analogous  plot  of  half-power  width.  In  both  figures,  a  and  L 
are  expressed  in  units  of  the  acoustic  wavelength,  and  direct  arrivals  are  omitted.  Because  the  scatter¬ 
ing  cross-section  formulas  contain  only  the  ratios  of  these  parameters  to  the  wavelength,  the  calcula¬ 
tions  are  independent  of  frequency.  They  encompass  a  wide  range  of  values  for  a  and  L.  For  a  24- 
wavelength  horizontal  array  in  this  configuration,  the  curves  predict  good  performance  (ASG  within 
~2  dB  of  reference  level)  when  o/\  <  2  and  L/X  >  120.  At  100  Hz,  these  values  are  about  the 
same  as  those  used  in  the  Gorda  Rise  example. 


2700  m 
3000  m 


212  km 


Fig.  16  —  Source-receiver  configuration  for  the 

calculations  of  Figs.  18  to  21.  The  source  is  at  5;  the 
horizontal  array  is  centered  at  C . 
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Normalized  rms  Roughness  a/ A 


?ig.  17  —  Array  signal  gain  as  a  function  of  rms  bottom  roughness  for  several 
alues  of  isotropic  bottom-height  correlation  length.  Lengths  are  normalized  by 
wavelength;  only  the  scattered  energy  is  included  in  the  computation. 


Fig.  18  —  Half-power  beamwidth  as  a  function  of  rms  bottom  roughness  for 
several  values  of  isotropic  bottom-height  correlation  length.  Lengths  are 
normalized  by  wavelength;  only  the  scattered  energy  is  included  in  the 
computation. 
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Finally,  Figs.  19  and  20  illustrate  cases  where  the  bottom  height  spectrum  is  anisotropic,  i.e., 
the  correlation  lengths  Lx  and  Ly  in  the  mean  bottom  plane  are  different  (see  Appendix  D).  In  this 
context  the  x -direction  is  upslope  in  the  bottom  plane,  and  the  y -direction  is  transverse  to  the  slope. 
Again,  direct  arrivals  are  omitted.  The  plots  show  that  the  scattering  is  strongly  dependent  on  the 
transverse  correlation  length  and  virtually  independent  of  the  slopewise  correlation  length.  Since  nei¬ 
ther  bottom  nor  array  is  yawed  in  these  calculations,  the  upslope  direction  is  perpendicular  to  the 
array.  The  results  show  that  although  scattering  in  the  vertical  may  change  the  level  of  sound  radia¬ 
tion  reaching  a  horizontal  array,  it  does  not  affect  the  detected  directionality. 


Transverse  Correlation  Length  Ly/A 

Fig.  19  —  Array  signal  gain  for  an  anisotropic  rough  bottom.  Lx  and  Ly  are  the 
bottom-height  correlation  lengths,  parallel  and  perpendicular,  respectively,  to  the 
propagation  direction.  Lengths  are  normalized  by  wavelength;  only  the  scattered 
energy  is  included  in  the  computation. 


5.  SUMMARY  AND  CONCLUSIONS 

We  have  described  a  practical  model  that  predicts  statistically  the  effect  of  rough-bottom  acous¬ 
tic  scattering  in  situations  where  one  bottom  interaction  is  dominant.  The  ORBS  model  permits  the 
use  of  measured  environmental  data  for  the  sound-speed  structure  and  bottom-roughness  statistics,  so 
that  realistic  calculations  can  be  made  for  specific  ocean  areas.  The  bottom  is  characterized  in  the 
model  by  a  basic  and  widely  used  measure  of  roughness,  the  power  spectral  density  of  residual 
bottom-height  variations  in  the  spatial  frequency  domain.  This  roughness  spectrum  requires  only  a 
few  parameters,  and  its  functional  form  in  the  program  can  readily  be  changed  to  conform  to  experi¬ 
mentally  measured  spectra.  Refraction  in  the  water  column,  bottom  reflection  and  scattering  in  three 
dimensions,  and  (optionally)  anisotropic  bottom  roughness  are  all  considered  in  the  model. 

The  report  also  demonstrates  the  application  of  an  exponential  substitution  method  that  provides 
an  accurate  and  computationally  efficient  approximation  for  the  Kirchhoff  scattering  integral  for  large 
values  of  rms  roughness.  The  approximation  takes  about  1%  of  the  time  required  for  numerical 
evaluation  of  the  integral.  With  this  approximation,  a  problem  involving  a  total  of  about  10-6  scatter¬ 
ing  paths  requires  about  30  min  of  CPU  time  on  the  VAX-11/780  computer.  (The  present  code  does 
not  use  an  array  processor.) 
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Fig.  20  —  Half-power  beamwidth  for  an  anisotropic  rough  bottom.  Lx  and  Ly 
are  the  bottom-height  correlation  lengths,  parallel  and  perpendicular, 
respectively,  to  the  propagation  direction.  Lengths  are  normalized  by 
wavelength;  only  the  scattered  energy  is  included  in  the  computation. 


The  short-range  sample  calculation  in  Section  4  shows  the  use  of  the  model  in  a  realistic  Pacific 
Ocean  environment  with  known  bottom  roughness  characteristics.  ORBS  computed  the  ensemble- 
averaged  directional  distribution  of  energy  that  is  forward-scattered  into  the  region  between  the  source 
and  the  first  strong  bottom  reflections.  For  a  smooth  reflecting  bottom,  very  little  energy  reaches  this 
region.  When  incoherent  scattering  is  included  in  the  calculation,  the  peak  angular  level  at  the 
receiver  has  dropped  only  —  8  dB  at  a  point  4  km  within  the  region.  The  ability  to  make  this  stochas¬ 
tic  calculation  in  three  dimensions  makes  ORBS  a  unique  predictive  model.  Average  directional  dis¬ 
tributions  can  be  calculated  with  deterministic  models  by  Monte  Carlo  methods,  but  this  requires  con¬ 
siderably  more  computer  effort,  especially  if  done  three-dimensionally. 

The  second  sample  case  shows  how  the  model  can  be  used  for  practical  array  design.  For  a 
particular  configuration  of  source,  ocean  bottom,  and  receiver,  the  results  of  several  ORBS  runs  are 
used  to  produce  graphs  relating  array  performance  to  the  roughness  parameters.  For  example,  for  a 
360-m  array  at  100  Hz,  array  signal  gain  is  within  2  dB  of  reference  level  when  the  rms  bottom 
roughness  a  <  30  m  and  the  correlation  length  L  >  2000  m.  Several  other  combinations  of  array 
and  bottom  roughness  parameters  can  be  used  for  such  studies  to  evaluate  the  effect  of  ocean  bottom 
roughness  on  array  performance. 

The  methods  used  in  this  three-dimensional  boundary-scattering  model  can  be  applied  to  other, 
more  general,  problems.  One  multiple-bounce  situation  is  underice  sound  propagation  in  the  Arctic 
where  the  rough  boundary  is  at  the  surface  and  the  bottom  is  not  encountered;  the  coherence  function 
approach  can  be  applied  here  also. 
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Appendix  A 

COHERENCE  FUNCTION  DEFINITIONS  AND  PROPAGATION 


In  this  appendix  we  derive  some  explicit  results  for  computing  the  transformed  coherence  func¬ 
tion  (or  “angular  power  distribution”)— how  it  propagates  in  the  absence  of  surface  interaction  and 
how  it  changes  upon  encountering  a  rough  interface.  Most  discussions  of  the  coherence  function  rely 
on  the  parabolic  approximation  to  the  Helmholtz  equation,  thereby  choosing  a  preferred  direction  of 
propagation.  Since  in  the  bottom-scattering  problem  there  are  two  “natural”  directions,  the  horizon¬ 
tal  propagation  direction  and  the  normal  to  the  mean  scattering  surface,  we  couch  our  discussion  in 
terms  of  the  Helmholtz  equation  to  avoid  picking  a  preferred  direction.  The  resulting  coherence  func¬ 
tion  can  then  be  computed  with  respect  to  arbitrary  planes. 

In  general,  the  two-point  mutual  coherence  function  T(r,  s)  is  given  in  terms  of  the  complex 
pressure  field  p  (r)  by 


T(r,  s)  =  <p( r  +  y)p*(r  -  y)>  .  (Al) 

The  angular  brackets  indicate  an  average  over  the  ensemble  of  ocean  boundaries,  sound-speed  pro¬ 
files,  or  sources.  To  begin,  we  assume  that  our  ensemble  is  degenerate,  having  only  one  member. 
Since  the  pressure  p  obeys  the  Helmholtz  equation 

V2  p(r  +  +  kl  n2( r  +  y)  p(r  +  y)  =  0 

V2p*(r  -  y)  +  &o  «2(r  -  y)/>*(r  -  y)  =  0,  (A2) 

we  find  an  equation  for  F  by  multiplying  the  first  of  these  equations  by  p*{ r  —  s/2)  and  the  second 
by  p(r  +  s/2)  and  then  subtracting.  The  result  is 


•  -f  r(r,  s) 
3s 


-*o 


n2( r  +  —)  -  n\ r  -  —) 


r  (r,  s) 


~  —  lefts  ■  Vn2  (r)T(r,  s).  (A3) 

The  last  approximation  is  called  the  “locally  quadratic”  approximation.  It  is  valid  when  «2(r)  is 
slowly  varying  over  the  domain  of  the  separation  variable  s  for  which  T  is  significant.  If  n2  is  a 
quadratic  function  of  r,  the  locally  quadratic  approximation  is  exact. 
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The  solution  of  Eq.  (A3)  is  found  by  Fourier  transforming  with  respect  to  the  separation  vector 
s  in  three  dimensions.  (This  differs  from  the  more  usual  procedure  of  transforming  s  over  a  plane.) 
We  then  find  that  in  the  locally  quadratic  approximation  the  transformed  coherence  function 

r  (r,  0)  =  j  e  ,k°e  s  r  (r,  s)  d3s  (A4) 

satisfies  the  equation 


0  '  ~t  f  (r’  6)  +  2  Vn2(r)  M  f  (r>  e)  =  °'  (A5) 

This  first-order  linear  equation  is  solved  by  the  method  of  characteristics.  In  this  case  the  characteris¬ 
tics  turn  out  to  be  the  rays  of  geometric  acoustics,  satisfying  the  equations 


30  1 

-  =  -  W(r). 


(A6) 


The  variable  t  is  related  to  the  arc  length  of  the  ray  by 

ds  =  n(r)dr.  (A7) 

From  the  ray  equations,  Eq.  (A6),  it  follows  that  “energy”  is  conserved  along  the  characteristics. 
That  is,  if  the  pairs  (r,0)  and  (r0,  0O)  are  joined  by  a  ray,  then 

e2  -  n2i r)  =  6q  -  nz( r0).  (A8) 

The  vector  0  is  tangent  to  the  rays  and  gives  the  local  direction  of  energy  propagation.  The  signifi¬ 
cance  of  the  characteristics  for  the  coherence  function  is  that  Eq.  (A5)  implies  that  f  is  constant  along 
characteristics,  just  as  is  the  “energy”  of  Eq.  (A8). 

It  is  tempting  to  interpret  T  as  the  joint  distribution  of  energy  in  space  and  in  direction  because 
the  marginals  of  T  do  indeed  give  the  spatial  and  angular  distribution  of  energy.  However,  since  T 
can  be  negative  for  some  pressure  fields,  such  an  interpretation  should  only  be  used  with  great  cau¬ 
tion.  What  is  more  certain  is  that  f  contains  the  acoustic  information  necessary  to  determine  quad¬ 
ratic  functionals  of  the  acoustic  field  such  as  outputs  of  a  square-law  detector. 

Equation  (A5)  implies  that  T  is  constant  along  the  rays  of  geometric  acoustics.  Hence 

f (r,  0)  =  f(r0,  0O)  (A9) 


where  (r,  0)  and  ( i'o ,  0g)  are  connected  by  a  ray.  Note  again  that  0  is  a  three-dimensional  vector; 
when  the  parabolic  approximation  is  used,  the  vector  corresponding  to  s  is  a  two-dimensional  vector 
lying  in  a  plane  normal  to  the  principal  direction  of  propagation.  Here  the  magnitude,  rather  than  the 
direction  of  0,  is  restricted  by  Eq.  (A8).  To  account  explicitly  for  this  restriction,  we  define 

7  (r,  0;  E)  =  8  (02  -  n2( r)  -  E)  f  (r,  0).  (A10) 
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Since  02  —  n2( r)  is  constant  along  rays  as  is  f ,  it  follows  that  7  is  constant  along  rays.  Furthermore. 
Eq.  (A  10)  implies  that  T  is  given  in  terms  of  7  by 

f (r,  0)  =  J  7  (r,  0;  E)  dE  (All) 


It  often  happens  that  information  about  the  pressure  field  and  hence  f  is  restricted  to  a  plane 
surface.  To  restrict  s  to  lie  in  a  particular  surface,  we  can  integrate 


Jv  0  r 

27  J  r  (r,  0)  dex  =  \e 


<p(r  +  y)p*(r  -  y)>  d20 p 


=  r(r,  0,) , 


(A  12) 


where  sx  is  the  component  of  s  perpendicular  to  the  surface  in  question  and  Sy  is  the  component 
parallel  to  that  surface.  Equations  (A  10)  and  (All)  can  be  used  to  write  the  restriction  of  T  to  a  sur¬ 
face  as 


A,  *0  r“  f(r0,«0|)|»oi  =  V^Tn2(r0)  -  9^|| 

r<r’#l|)  ■  2?  J-  T  nH r)  -  9. 


By  using 
this  equation  becomes 


=  ko_  f(r.  0|,)|0X  =  Vk  +n2(r)  -  0f 
27r  _0°  2V£  +  n2(r)  -  0j 

E  =  0O2±  +  -  «2(r0)  and  </£  =  20Ol  d60x  , 

t,,  c  f\r0,  0O|  1 0o_L)0o_Lfi?0o 

F(r’  *l|)  =  TT  J  ^  -  T,  ~~T 


2ir  v  ^oil  +  -  «2(r0)  +  n2(r)  -  ' 


(A  13) 

(A14) 

(A  15) 


To  the  extent  that  f(r0,  0O)  is  sharply  peaked  about  0q||  +  ^o±  —  n 2(r0)  =  E,  we  can  write  Eq. 
(A  15)  as 


Vf  -  dL  +  n 2(r0)  „ 

r(r,  0||)  =  7= — a - =£-  r(r0,  0O|). 

V£  -  0?  +  n2(r) 


(A  16) 


When  f  is  developed  from  the  parabolic  approximation,  both  0  and  0O  are  associated  with  normals  to 
vertical  planes.  In  the  present  development,  the  plane  associated  with  0  need  not  be  parallel  to  that 
associated  with  0O,  and  neither  is  required  to  be  vertical.  In  practice  E  is  taken  to  be  zero. 

Scattering  of  the  Coherence  Function 

When  a  rough  surface  is  encountered,  the  angular  distribution  of  energy  is  suddenly  altered  and 
the  distribution  of  scattered  energy  is  a  functional  of  the  distribution  of  incident  energy.  We  therefore 
make  the  assumption  that  the  coherence  function  transform  associated  with  the  scattered  energy  is 
given  by 


q-  «  rscat  (rfi,  q)  =  j  £  (k  -  q)k-  n  rinc  (rB,  k)  k03  d?  k.  (A17) 
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The  transition  kernel  £  is  closely  related  to  the  differential  cross  section  per  unit  area.  This  descrip¬ 
tion  of  surface  scattering  is  local  in  the  sense  that  the  scattered  f  at  the  point  FB  on  the  mean  scatter¬ 
ing  plane  is  determined  by  the  incident  F  at  that  same  point.  In  general  this  is  an  approximation 
whose  validity  is  hard  to  assess.  We  also  note  that  the  definition  of  T  is  problematic  when  the  propa¬ 
gation  region  is  bounded,  since  the  Fourier  transform  that  defines  f  extends  over  all  space  while  the 
pressure  field  is  defined  only  in  a  limited  region  of  space.  One  therefore  must  argue  that  f  is 
nonzero  only  over  a  limited  range  of  the  separation  variable  s,  or  that  only  a  limited  range  of  s  affects 
the  transform.  Assuming  the  vectors  r  +  s/2  lie  in  the  propagation  region,  the  infinite  transform  of  s 
presents  no  real  problems.  This  means  that  r  must  be  kept  away  from  the  boundaries. 

In  Eq.  (A17),  the  scattering  kernel  W  can  be  found  by  considering  what  happens  to  a  plane 
wave.  We  assume  that  F  can  be  defined  in  a  region  near  enough  to  the  scattering  surface  so  that  in 
this  region  the  sound  speed  may  be  considered  to  be  a  constant.  We  need  a  region  that  is  large 
enough  so  that  there  is  no  difficulty  in  the  definition  of  T  but  small  enough  that  refraction  can  be 
ignored.  If  such  a  region  exists,  we  consider  the  incident  plane  wave 

P(  r)  =  (A18) 


for  which  the  corresponding  function  f  is 

3 

S3  (0  -  k)  (A  19) 

=  Tine  (r,  0). 

The  corresponding  scattered  field  above  the  highest  points  of  the  scattering  surface  can  be  written  as  a 
superposition  of  plane  waves 


fW  (r>  0)  — 


2tt 

ko 


•^scat  (r)  J  e  q  eq±  1  A  (qy,  ky)  d2q y ,  (A20) 

where  q  ±  =  Vkfi  ~qf  ,  and  the  subscripts  ||  and  ±  refer  to  components  parallel  to  and  perpendic¬ 
ular  to  the  mean  scattering  surface.  If  only  a  finite  portion  of  the  scattering  surface  is  rough,  the 
scattering  amplitude  in  the  expression 

J’scat  (r)  =  f  (?,  k;  k0)  (A21) 

is  related  to  the  plane  wave  amplitudes  /4(k,  q)  by 

/(q,  k;  k0)  =  2tt k0  q  x  4(qy,  ky).  (A22) 

The  differential  cross  section  a  for  such  a  finite  scattering  patch  is  given  by 

a(k  -q)  =  -LAaJSlL  _  (2tt)2  q2±  |A(qN,  ky)|2  (A23) 

K0 
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from  Ref.  Al.  If  p  such  roughness  elements  per  unit  area  scatter  independently  and  incoherently 
then  pa  is  the  differential  cross  section  per  unit  area.  Correspondingly,  if  the  surface  is  randomly 
rough  with  homogeneous  statistics,  the  fluctuations  of  the  plane-wave  scattering  amplitudes  satisfy 

<A A  (0,  +  -^kk,,)  AA*  (0„  -  -I1,  k|,)>  =  52  (q„)  A  (0, ,,  k,,).  (A24) 


The  function  f1  that  arises  from  products  of  the  amplitude  fluctuations  A A  ,  rs‘"atoh,  is  given  by 

j  g‘>rs«  ei4,±s±)  A  (0|,  k|j)  d<t>\\  ds  ,  (A25) 


rsr  (r,  0)  =  J  e‘ 


where  (f>±  =  —  <t>j,  and  where  we  assume  that  |k|  |  <  k0.  It  follows  from  Eq.  (A25)  that 

f1  scat°h  is  given  by 

r‘cnat0h  (r,  e)  =  (27T)3  5(1  ~  |g|)-  Vl  -  e\  A  (*0  «|,  k«)-  (A26) 

Kq 


On  the  other  hand,  if  in  Eq.  (A17)  we  use 


Fine  (r.  4>) 


53  (0  -  k), 


it  follows  that 


0  ■  n  fs'"“h  (r,  0)  =  (27 r)3  (k  ft)  £incoh(k  -  0).  (A27) 

Hence  we  conclude  that  the  scattering  kernel  for  incoherent  scattering  is  given  in  terms  of  fluctuations 
of  the  scattered  plane  wave  amplitudes  by 


q-  n  ^  (k  -  k0  0)  =  e\  5(1  ~  A  (k0  k„) .  (A28) 

Kq 

By  comparing  to  results  for  a  surface  consisting  of  distinct,  incoherent  scattering  patches,  we  find  that 
A,  the  strength  of  the  plane-wave  amplitude  fluctuations,  is  related  to  the  differential  cross  section  per 
unit  area  by 

A  (p n ,  <?ii)  =  a  (q  -  p).  (A29) 

PI 

There  is  also  a  part  of  f  associated  with  the  average  of  the  plane  wave  scattering  amplitudes, 
<A( q,  k)>.  For  homogeneous  surface  statistics,  <A(q,  k)>  is  proportional  to  a  delta  function 

<A(q,  k)>  =  6(k||  -  q|)  a  (^||).  (A30) 

In  the  Kirchhoff  approximation,  for  example,  a  (qA  is  given  by 

«(q||)  =  e~2a2ql  .  (A31) 
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The  contribution  to  the  scattered  function  f  arising  from  the  coherent  scattered  amplitudes  is 

rS5(r,«)  =  |a(0)|2finc(r,  0inc) 


or 


rscc°at  (r,  ®|,  ex)  =  |a(0)|2  finc  (r,  0„,  -0X).  (A32) 

Energy  conservation  requires  that  the  incident  energy  flux  equal  the  total  of  the  scattered  energy  flux 
normal  to  the  mean  scattering  plane.  In  terms  of  a  and  A  this  means 

£ 

lfl(?)l2  + 1  ^7  A  ^n*  qiP  d2k\\  =  l-  (A33> 


In  the  Kirchhoff  approximation,  this  relationship  is  only  approximately  satisfied. 

Reference 

Al.  D.H.  Berman  and  J.S.  Perkins,  “The  Kirchhoff  Approximation  and  First-Order  Perturbation 
Theory  for  Rough-Surface  Scattering,”  J.  Acoust.  Soc.  Am.  78,  1045-1051  (1985). 
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Appendix  B 

ARRAY  SIMULATION  AND  THE  COHERENCE  FUNCTION 


In  this  appendix,  we  show  how  array  response  is  expressed  in  terms  of  the  three-dimensional 
coherence  function  transform  discussed  in  Appendix  A.  We  confine  our  discussion  to  a  linear  array, 
but  the  array  may  be  arbitrarily  oriented.  The  possibility  of  off-broadside  orientation  requires  use  of 
the  three-dimensional  coherence  function  transform. 

Let  c  be  a  three-vector  that  locates  the  center  of  a  linear  array  of  hydrophones.  Let  eL  denote  a 
unit  vector  pointed  along  the  array.  Then,  the  signal  that  results  from  (a)  steering  the  array,  (b)  shad¬ 
ing  the  array,  and  (c)  summing  the  steered,  shaded  fields  at  each  hydrophone  is  given  by 

A(6S)  =  j  W(y)  p( c  +  ye/>"'V  sine‘dy  .  (Bl) 

The  shading  function  W  allows  for  either  discrete  or  continuous  arrays.  The  pressure  p  at  the  point 
c  +  yeL  is  denoted  as  p(c  +  eLy).  Often  the  signal  out  of  an  array  is  squared  to  form  the  power. 
Equation  (Bl)  implies  that  the  power  is  given  by 

|A(0J|2  =  n  ^Or)  W*(y2)p(c  +  y}eL)P*(c  +  yi  h)e  ^  >2>  sm  $s  dyxdy2.  (B2) 

This  expression  contains  the  product  of  pressures  at  two  different  points  and  therefore  can  be 
expressed  in  terms  of  f .  Changing  to  mean  and  difference  coordinates,  we  can  write 


P(c  +  yieL)p*(c  +  y2eL) 


2 

f  (C  + 


y\ 


eiko(y~y2)h 


Ldye. 


(B3) 


Therefore  the  array  power  A 2  is  given  by 


|  A(0S)  |2  =  ds  W(y  +  ±)W*(y-  j) 


2ir 


ld36f(c+yeL,e)e 


fays  (6  ■  eL  —  sin  $s ) 


(B4) 


A  more  succinct  expression  for  the  power  can  be  obtained  by  defining  an  analog  of  f  for  the  shading 
functions.  Let  a (y ,  0,  6S)  be  given  by 


a(y,6,6s)  =  \  W(y 


+  y)  W*(y  -  y)  e 


‘V(®£ i-sinfls)  ^ 


(B5) 


Equation  (B4)  becomes 


l^(^)|2  =  H  fl(y,0,0i)f(c  +  ytL,  0)d3e  dy  . 


(B6) 
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Note  that  a  (Y ,0,0,)  depends  on  0  only  through  eL  •  0  =  sin  The  angle  $  is  the  arrival  angle  at  the 
array  measured  from  broadside.  It  follows  that  Eq.  (B6)  simplifies  to 

\A(6s)\2  =  |  j  a(y,\J/,es)f( c  +  yeL,6)d3d  dy  .  (B7) 


We  assume  that  by  judicious  choice  of  k0  we  can  make  f  sharply  spiked  about  |0|  =1.  By  using 
Eqs.  (A10)  and  (All),  we  find 

\A@s)  1 2  =  j$a(y,t,0s)  ^-f(c  +yei,\0,n,0  =  N/E  +^)cos\[/d\pdQdy  dE  ,  (B8) 

where 


n 2  =  n\  c  +  yeL). 

If,  in  fact,  f  is  proportional  to  a  delta  function  in  E,  f  =  f0  5  (E),  then 

\A(6s)\2  =  H  f  -^n2(c  +  yeL)a(y,\P,ds)r 0(c  +  yeL,\0,fl)  cos  dy  (B9) 

where 


ro(r,\0,Q)  =  f0  (r,n(r)  (sin  eL  +  cos  \p  cos  fi  e//  -f-  cos  \j/  sin  Q  e^)).  (BIO) 

(The  unit  vectors  e#  and  i?y  are  defined  below.) 

We  model  a  shaded  array  of  N  discrete  hydrophones  located  at  positions  along  the  array,  yf,  by 

W(y)  =  £  Wt  8  (y  -  y;).  (Bll) 


The  shading  function  applied  to  each  hydrophone  is  denoted  wt.  The  function  a(y  ,^,0f)  that  sum¬ 
marizes  intrinsic  array  properties  is  now  given  by 

=  I ;s(y  -  -  *  yj )  (B12) 

ij  2 

Note  that  a  is  nonzero  not  only  at  hydrophone  positions  but  also  between  hydrophone  positions  at 
CV;  +  yj)/ 2.  By  using  this  expression  for  a ,  the  array  power  can  be  written 

\M0S)\2  = 


Elly  ”2(c  +  ?  -  ®L)fl(y>'/'>0,)r'o(c  + 

U  Z 


T/  +  y< 

- - — -eL,\0,Q)  cos  \p  d\p  d fi. 


(B13) 


Equation  (B13)  is  the  central  result  of  this  appendix. 
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The  function  a  summarizes  the  intrinsic  array  properties.  We  now  comment  on  the  description 
of  the  array  geometry.  We  have  denoted  a  unit  vector  pointing  along  the  array  by  eL.  Let  z  be  the 
unit  vector  pointing  downward.  We  define  a  unit  horizontal  vector  perpendicular  to  the  array  by 


eL  X  z 
\h  x  z| 


(B14) 


If  there  is  no  array  yaw,  eH  points  generally  in  the  direction  from  the  receiver  back  to  the  source.  A 
third  orthogonal  array  fixed  coordinate  is  defined  by  eH  x  eL  =  ev.  The  angles  0  and  12  are  chosen 
so  that  an  arbitrary  direction  of  arriving  energy  0  can  be  written 

0  =  sin  0  eL  +  cos  0  cos  12  eH  4-  cos  0  sin  fi  .  (B15) 


If  the  array  is  yawed  by  an  amount  a  and  then  tilted  by  an  amount  /3  (order  is  important  for 
rotations!),  then  eL  is  given  by 


eL  =  —  sin  a  cos  /?  x  —  cos  a  cos  /3y  +  sin  /3  z. 


(B16) 


In  similar  fashion,  the  array-fixed  unit  vectors  tH  and  eK  can  be  expressed  in  terms  of  a,  (3,  and  the 
space-fixed  unit  vectors  x,  y,  z. 

If  energy  arrives  along  the  unit  vector 

0  =  cos  0  cos  <j>  x  —  sin  0  cos  <j>  y  +  sin  <j>  z ,  (B17) 

then 

eL  ■  0  =  cos  0  cos  0  sin  (0  —  a)  +  sin  0  sin  0  (B18) 


=  sin  0'. 


Thus  if  a  =  0,  sin  0  =  •  0  =  sin  0'.  This  is  Eq.  (11)  of  the  text.  Similar  considerations  are 

required  to  map  directions  0  and  12  specified  relative  to  array-fixed  axes  into  angles  relative  to  space- 
fixed  axes.  Finally  note  that  the  angle  defined  in  Eq.  (B15)  is  the  cone  angle  appearing  in  Eq.  (B13) 
and  in  Fig.  2. 


* 
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Appendix  C 

TRACING  THE  CHARACTERISTIC  PATHS 


The  tracing  algorithm  used  in  ORBS  computes  closed-form  solutions  for  layer-to-layer  path 
increments  in  a  range-independent  stratified  sound-speed  structure,  so  that  the  paths  may  be  computed 
efficiently  for  a  general  multilayer  deep-ocean  profile.  The  profile  is  entered  as  a  set  of  depth-sound 
speed  pairs  (z,- ,  c, } ,  and  an  interpolating  function  of  the  form 

c(z)  =  c{[  1  -  2 gt(z  -  z,)]“1/2,  Z;  <  z  <  Z;+i  (Cl) 

is  used;  g,  is  a  fitting  parameter  for  the  ith  layer.  One  way  of  expressing  the  solution  for  the  path 
for  this  interpolating  function  [Cl]  is 


tan  <j>  =  tan  </>,  -  rgt  (cm/c,)2,  (C2) 

where  <t>  is  the  angle  measured  down  from  the  horizontal,  r  is  the  range  in  the  ith  layer  (the  path 
enters  the  layer  with  angle  <£,),  and  the  parameter  cm  is  defined  by  Snell’s  Law: 

c  =  cm  cos  <t>.  (C3) 


In  addition  to  the  layer-by-layer  trace  of  path  increments,  we  need  an  efficient  way  of  tracing  a 
path  to  its  intersection  with  the  bottom.  In  ORBS,  this  is  done  by  assuming  that  the  paths  contact  the 
bathymetry  only  in  deep  water  ( z  >  3000  m),  where  the  depth  gradient  of  the  sound  speed  is  approx¬ 
imately  constant  and  a  simplified  tracing  algorithm  can  be  used.  The  bottom  layer  of  the  profile  is 
fitted  to  the  linear  interpolating  function 


c  (z )  —  C  +  Gz  . 


(C4) 


Then,  the  path  solution  can  be  written 


sin  <t>  =  sin  0,-  -  rG/cm  ,  (C5) 

where  0,  </>,-  and  cm  are  as  defined  above,  and  r  is  the  range  in  the  bottom  layer.  Now  suppose  the 
mean  bottom  plane  intersects  the  vertical  plane  of  the  path  in  a  line  whose  equation  is 

z  =  z0  +  hr .  (C6) 

To  find  the  point  where  the  path  intersects  the  bottom  plane,  define  cb  as  the  sound  speed  at  the  bot¬ 
tom  intersection  and  combine  the  last  three  equations  (and  Snell’s  Law)  to  obtain  a  quadratic  in  cb : 

(/  +  h2)c£  -  2 Bcb  +  B2  —  (l hcj 2  =  0,  (C7) 
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where  B  =  C  +  Gz0  +  hcm  sin  </>,■ .  This  equation  can  be  solved  by  the  quadratic  formula.  (If  no 
solution  exists,  the  path  misses  the  bottom  and  is  not  a  contributor  to  the  scattered  or  reflected 
energy.)  The  range  and  depth  of  the  bottom  hit  are  now  linear  functions  of  cb . 

Reference 

Cl.  M.A.  Pedersen  and  D.F.  Gordon,  “Normal-Mode  Theory  Applied  to  Short-Range  Propagation 
in  an  Underwater  Acoustic  Surface  Duct,”  J.  Acoust.  Soc.  Am.  37,  105-118  (1965). 
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Appendix  D 

COMPUTING  THE  SCATTERING  CROSS  SECTION 


The  scattering  cross  section  governs  the  fraction  of  energy  impinging  upon  the  bottom  that  is 
scattered  in  each  direction.  It  is  a  function  of  the  bathymetric  spectrum  and  the  direction  of  the 
incoming  energy.  The  model  incorporated  in  ORBS  is  based  on  a  combination  of  small  perturbation 
theory  and  the  Kirchhoff  approximation  [Dl].  This  model  reduces  to  perturbation  theory  for  small 
surface  heights  and  to  conventional  Kirchhoff  theory  for  gently  undulating  surfaces  with  large  surface 
heights.  Thus,  the  region  of  validity  of  the  model  is  at  least  as  great  as  the  union  of  the  regions  of 
validity  of  its  limiting  forms. 

Let  q  be  the  wave  vector  incident  upon  the  bottom  and  let  k  be  the  wave  vector  in  a  scattered 
direction.  According  to  Kirchhoff  theory,  the  proportion  of  the  incident  energy  reflected  coherently 
(in  the  specular  direction)  is  given  by 


^refl 


e 


“4  fcV 


(Dl) 


where  qz  is  the  component  of  q  normal  to  the  mean  surface,  and  a  is  the  rms  residual  bottom  height. 
The  remainder  of  the  incident  energy  is  scattered  in  the  nonspecular  directions.  The  cross  section  per 
unit  area  for  incoherent  scattering  from  the  incident  direction  q  to  the  scattered  direction  k  may  be 
written  [Dl]: 


a(q  -  k) 


— X—  7M  +  2R  M-l'l.-l  +  b-l  7-l,-l) > 

16ir  qz 


(D2) 


where  R  is  an  elastic  coefficient  between  —  1  (perfectly  soft  bottom)  and  + 1  (perfectly  hard  bottom) 
and 

=  J  exp  [i (q  -  k)  r]  (exp  [Xm cr2 kK(/- )]  -  1  }/{\m\n)d2r  ,  (D3a) 

Xm  =  q  +  mk,  m  =  ±  1 ,  (D3b) 

and 

brn  =  [(4z  +  mK )2  +  I  q||  —  k||  |2]  exp  (-  X2  a2 /2) .  (D3c) 

In  the  Eqs.  (D3),  the  subscript  z  designates  the  component  of  a  vector  normal  to  the  mean  surface, 
and  the  subscript  ||  designates  the  components  parallel  to  the  surface.  The  function  W(r)  in  Eq. 
(D3a)  is  the  normalized  correlation  function  of  surface-height  fluctuations, 

W(r)  =  <h (r )h (0) > / < h 2(0) >  ,  (D4) 
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where  h  (r)  is  the  height  of  the  scattering  surface  above  the  point  specified  by  the  two-dimensional 
vector  r  in  some  reference  plane.  The  two-dimensional  surface  height  spectrum  is  defined  by  a  two- 
dimensional  Fourier  transform  of  W{ r): 

5(k)  =  a2  |  j  exp  (ik-  r)  W(r)d2r  ,  (D5) 

where  the  rms  surface  roughness  <h\ 0)>  is  denoted  by  o2.  By  inversion  of  the  transform  of  Ea. 
(D5),  IT(r)  satisfies 

W(r)  =  (2ira)~2  exp  (-ik-  r)  S(k)  d2k  .  (D6) 


In  the  examples  to  follow,  consider  the  isotropic  bathymetric  correlation  function 

W(r)  =  ( r/Lc)Kx(r/Lc ),  (D7) 

where  Kx  is  the  first-order  modified  Bessel  function  of  the  second  kind,  as  a  function  of  the  scalar 
separation  r  =  |  r  | .  This  choice  of  W{r)  corresponds  to  a  k~4  dependence  in  the  high-wave-number 
region  of  the  spectrum;  specifically,  the  spectrum  obtained  by  transforming  Eq.  (D7)  is 

S(k)  =  4tt(oLc)2  [1  +  (kL )2]  ~2 ,  (D8) 

where  k  =  |  k  |  and  Lc  is  the  correlation  length  of  the  bottom  height  variations  about  the  reference 
plane.  (Other  choices  for  W{r)  may  be  considered,  such  as  a  Gaussian,  which  makes  S(k)  a  Gaus¬ 
sian  also.) 

Figure  D1  shows  the  results  of  cross-section  calculations  with  the  correlation  function  and  spec¬ 
trum  of  Eqs.  (D7)  and  (D8),  with  Lc  =  1600  m,  a  =  30  m,  and  for  a  100  Hz  source  impinging  at  a 
grazing  angle  of  20°  to  the  mean  bottom.  For  such  a  large  roughness,  we  would  not  expect  perturba¬ 
tion  theory  to  be  accurate  but  would  expect  Kirchhoff  theory  to  be  valid.  The  cross  section  is  plotted 
in  Fig.  Dl,  where  the  modified  Kirchhoff  result  is  plotted  on  the  top  as  a  function  of  outgoing  verti¬ 
cal  and  azimuthal  angles,  and  the  corresponding  perturbation  and  standard  Kirchhoff  results  are 
below.  The  plot  also  shows  the  position  and  value  of  the  maximum  of  each  of  the  three  formulations. 
As  expected,  the  modified  Kirchhoff  calculation  is  accurate  in  this  region  (as  shown  by  comparison 
with  the  standard  Kirchhoff  result),  while  perturbation  theory  is  in  error  by  several  orders  of  magni¬ 
tude.  Further  examples  are  contained  in  Ref.  D2. 

As  stated  in  Section  2,  hundreds  of  thousands  of  evaluations  of  the  cross  section  are  needed  in  a 
typical  ORBS  calculation.  The  most  time-consuming  part  of  the  calculation  is  the  numerical  computa¬ 
tion  of  the  integral  in  Eq.  (D3a),  which  is  needed  for  Eq.  (D2).  We  have  simplified  the  calculation 
by  a  technique  that  falls  into  the  class  of  approximations  called  exponential  substitution  [D3,D4],  The 
integral  to  be  evaluated  is  of  the  form 

7(Q)  =  j  J  exp  ( iQ  •  s)  [  exp  [HW(s)]  -  1  }d2s  .  (D9) 

We  assume  that  we  can  approximate  this  integral  by 

J(Q)  =  j  J  exp  (iQ  -  s)  \W(s/s0)d2s  ,  (D10) 

where  A  and  s0  are  parameters  to  be  chosen.  With  a  change  of  variables  we  obtain 

•7(Q)  =  So  A  j  j  exp  (isQ- x)  W(x)d2x  .  (Dll) 
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0 

* 


MODIFIED  KIRCHHOFF 
S  (  20.0,  0.0)  =  43.08 


100  Hz 

=  20  deg 
a  =0  deg 

a  =  30  m 
L,  =  1600  m 
l_2  =  1600  m 

HARD  BOTTOM 


PERTURBATION  RESULT 
£  (  20.0,  0.0)  =  999999.00 


KIRCHHOFF  RESULT 
£  (  20.0,  0.0)  =  43  08 


Fig.  D1  —  Scattering  cross  section  as  a  function  of  vertical  and  azimuthal 
angles  obtained  by  using  the  exact  formulation 
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The  integral  in  Eq.  (Dll)  is  the  two-dimensional  Fourier  transform  of  the  bathymetric  correlation 
function,  which  by  definition  gives  the  bottom-height  spectrum  S  (Eq.  (D5)),  that  is, 

J (Q)  =  A  S(Q^0)/o2  (D12) 


We  now  choose  A  so  that  the  integrands  of  /  and  J  are  equal  for  s  =  0,  obtaining 

A  =  exp  (H)  -  1 .  (D13) 


There  are  many  ways  to  choose  s0.  It  is  convenient  and  accurate  to  require  that 

J  I  <2=0  =  J  I  e=o- 


(D14) 


This  means  that  the  approximation  is  exact  in  the  specular  direction.  The  result  is 


S  n  — 


a2  |  |  exp  [tfJT(s)  -  1  ]d2s 


sme 


H 


1) 


(D15) 


This  defines  j0  as  a  function  of  H  that  in  general  is  very  smooth  and  readily  tabulated.  Its 
evaluation  is  then  simply  a  matter  of  table  look-ups  and  interpolation.  Because  as  S(k )  is  often  given 
in  a  simple  closed  form  (for  example,  if  W  is  Gaussian  or  as  in  Eq.  (D7)),  the  evaluation  of  the  cross 
section,  Eq.  (D2),  is  greatly  reduced  in  complexity  by  exponential  substitution. 

Figure  D2  shows  the  result  of  using  the  approximation  for  the  same  case  as  Fig.  Dl.  The 
differences  in  the  modified  Kirchhoff  result  between  Figs.  Dl  and  D2  are  minimal.  The  peak  values 
(for  Q  =  q  —  k  =  0)  are  virtually  identical,  but  of  course  this  is  one  of  the  conditions  set  by  Eq. 
(D13).  The  difference  in  total  scattered  energy  corresponding  to  the  two  functions,  computed  by 
numerical  integration,  is  <7%.  On  a  VAX-11/780  computer,  it  took  90  min  of  CPU  time  to  calcu¬ 
late  the  results  in  Fig.  Dl,  while  Fig.  D2  took  —90  s,  with  a  significant  percentage  taken  up  in  the 
plotting.  Thus,  using  this  approximation  reduces  computational  time  by  almost  two  orders  of  magni¬ 
tude. 


The  bathymetric  statistics  are  easily  generalized  to  model  an  anisotropic  bottom.  An  example  of 
an  anisotropic  bottom  spectrum  with  different  correlation  lengths  in  the  x  -  and  y  -directions  is 

S(k)  =  4t r^L2  [1  +  k2L2]~2,  (D16) 

where  the  vectors  L  =  (Lx,  Ly)  and  k  =  (kx,  ky)  and  L2  =  LxLy .  This  spectrum  is  a  generalization 
of  the  spectrum  of  Eq.  (D7).  It  has  the  same  k~4  dependence  in  the  high-wave-number  tail  when 
evaluated  along  either  principal  axis  in  k -space.  By  use  of  the  coordinate  transformation 
k '  =  (kx  \fLx/Ly,  ky\lLy/Lx),  the  transform  integral  for  the  correlation  function  (Eq.  (D6))  can  put 
into  the  same  form  as  the  isotropic  case,  and  the  correlation  function  IT(r)  is  equal  to 

W(r')  =  (r'/Lxy)Kl(r'/Lxy),  (D17) 
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MODIFIED  KIRCHHOFF 
E  (  20.0,  0  0)  =  43.05 


100  Hz 

Pin  =  20  de9 
a  =  0  deg 


a  =  30  m 
L,  =  1600  m 
L2  =  1600  m 

HARD  BOTTOM 


pfrturbation  result 


KIRCHHOFF  RESULT 


E  (  20.0,  0  0)  =  999999.00 


Fig.  D2  —  Scattering  cross  section  as  a  function  of  vertical  and  azimuthal 
angles  obtained  by  using  the  exponential  substitution  approximation 
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where  the  vector  r  (x  / Lx ,  y  V Lx  / Ly ,  and  r'  —  |r'  | .  Then  the  Kirchhoff  integral  assumes 
the  same  form  as  in  the  isotropic  case,  and  the  same  method  of  exponential  substitution  outlined 
above  can  be  used  to  approximate  it.  The  result  is  that  J(Q)  is  equal  to 

J( Q)  =Ai02  SisoQ')/^,  (D18) 

where  the  vector  Q'  =  (Qx  \lLx /Ly ,  Qy  \lLy/Lx).  Anisotropic  bottom  statistics  are  thus  introduced 
by  specifying  two  correlation  lengths  in  the  principal  directions  and  using  them  to  define  a  scale 
change  in  the  coordinate  system  in  k -space.  This  type  of  anisotropic  spectrum  could  be  used  to 
represent,  for  example,  a  sloping  bottom  with  furrows  running  along  the  direction  of  the  slope.  This 
would  be  the  direction  with  larger  correlation  length. 
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Appendix  E 

PROGRAM  DATA  INPUTS 


This  appendix  consists  of  an  outline  of  the  input  data  for  ORBS  and  a  sample  data  file. 

There  are  two  ways  of  specifying  a  field  of  coherence-function  values  to  start  up  the  program: 

•  The  program  can  read  from  a  binary  file  of  values  of  the  coherence  function  transform,  com¬ 
puted  from  a  PE  computation,  for  a  PE  computation  for  a  succession  of  source-cylinder  ranges.  This 
allows  the  source-receiver  range  to  be  varied  for  a  given  computed  sound  field. 

•  A  simulated  field  can  be  used.  The  energy  at  the  source  cylinder  is  distributed  in  a  two- 
dimensional  Gaussian  about  a  mean  depth  and  a  mean  vertical  direction.  In  this  case  the  random- 
access  file  and  FORTRAN  unit  4  are  not  needed.  The  flag  for  choosing  the  startup  method  is  on 
record  4. 

Read  in  the  data  in  the  following  series  of  card  images  (note  that  it  was  considered  convenient  to 
use  the  format  (10F8)  for  all  data): 

Record  1:  Any  alphanumeric  title  in  columns  1-72 

Record  2:  RPKM,RCKM,ZC,ALPHAD,BETAD,XL,XN,XMODE 


RPKM 

—  Source  cylinder  range  (km) 

RCKM 

—  Range  (km)  of  center  of  array  (point  C) 

ZC 

—  Depth  (m)  of  point  C 

ALPHAD 

—  Array  yaw  (rotation  about  vertical)  (deg).  Zero 

yaw  means  the  array  is  broadside  to  the  source.  Positive 
yaw  is  counterclockwise  as  seen  from  above. 

BETAD 

—  Array  tilt  (deg)  in  vertical  plane  about  point  C. 

Positive  tilt  is  clockwise  as  seen  from  the  source. 

XL 

—  Length  of  array  (m) 

XN 

—  Number  of  phones  in  array 

XMODE 

—  A  flag  0.  to  7.:  set  4-bit  to  omit  direct  paths 

set  2-bit  to  omit  reflected  paths 
set  1-bit  to  omit  scattered  paths 
(But  enter  with  decimal  point) 
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Record  3: 


Record  4: 


SIGMAS, CLX,CLY,ZBC,ZMAX, SLOPED, YAWD,SPOW,GAMMAD 


SIGMAS 

CLX 

CLY 


ZBC 


—  rms  roughness  (m)  of  the  bottom;  i.e.,  the  standard 
deviation  of  residual  bottom  heights  about  mean  bottom  plane 

—  Correlation  length  (m)  of  residual  bottom  height 
distribution,  measured  in  slope  direction  in  mean  bottom  plane 

—  Correlation  length  (m)  of  bottom  height  distribution, 
measured  in  cross-slope  direction  in  mean  bottom  plane.  Note 
that  these  measurement  directions  can  be  rotated  by  specifying 

a  nonzero  skew  angle.  For  an  isotropic  bottom,  enter  CLY  equal 
to  CLX  (see  GAMMAD). 

—  Depth  of  mean  bottom  plane  (m)  directly  under  the  array 
center 


ZMAX  —  Maximum  water  depth  (m) 

SLOPED  —  Slope  (deg)  of  mean  bottom  plane.  A  positive 
slope  is  downward  from  receiver  toward  source. 

YAWD  —  Yaw  (rotation  about  vertical)  (deg)  of  bottom. 

Zero  yaw  means  isobaths  in  the  mean  bottom  plane  are 
broadside  to  the  source.  Positive  yaw  is  counterclockwise 
as  seen  from  above. 


SPOW  —  Inverse  power-law  exponent  for  bottom  roughness  spectrum. 

ORBS  models  a  power  spectral  density  for  residual  bottom 
heights  asymptotic  to  k**  (-SPOW),  where  k  is  the  spatial  wave 
number.  Values  of  4,  5,  and  6  may  be  used  for  power-law  spectra; 
a  value  of  0  will  cause  a  Gaussian  spectrum  to  be  used. 

(These  are  the  only  values  for  which  tables  for  the  exponential  substitution 
algorithm  have  been  provided.) 

GAMMAD  —  Roughness  ellipse  skew  angle  (deg)  (CCW  from  above). 

If  CLX  =  CLY,  the  value  of  GAMMAD  is  irrelevant  (enter  zero  in 
this  case). 

FHZ,  GAMO,  ZS,  SIGZS,  SBEAMD,  SIGSBD,  XNSVP 


FHZ  —  Acoustic  frequency  (Hz). 

GAMO  —  Peak  value  for  simulated  coherence  function  at  the  source 

cylinder,  as  described  above.  (Recommended  value:  1 .0  to  use 
simulated  source;  put  GAMO  =  0.  to  use  the  PE-generated  startup  file.) 

ZS  —  Depth  (m)  of  peak  coherence  function  value  for  simulated 

beamed  source 

SIGZS  —  Standard  deviation  in  depth  (m)  of  Gaussian  distribution 
for  simulated  beamed  source 
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SBEAMD  - 

Angle  (deg)  of  peak  coherence  function  value  for  simulated 

SIGSBD  — 

beamed  source.  AO  is  positive  for  a  beam  directed  downward  from  the 
source. 

Standard  deviation  in  angle  (deg)  of  Gaussian 
distribution  for  simulated  beamed  source 

(The  values  of  ZS,  SIGZS,  SBEAMD,  and  SIGSBD  are  irrelevant  if  GAMO  =  0.) 

XNSVP  —  Number  of  depth-sound  speed  pairs  in  profile  to  follow 
Record  5:  PSI1D,DPSI,PSI2D,  OMEGAlD,DOMEGAD,OMEGA2D,  DGRZD 


PSI1D  — 

Minimum  beam  angle  (deg)  for  calculation 

DPSID  — 

Beam  angle  increment  (deg) 

PS12D  — 

Maximum  beam  angle  (deg)  for  calculation 

OMEGA  ID  — 

Minimum  depression  angle  in  beam  cone  (deg) 

DOMEGAD  - 

Depression  angle  increment  (deg) 

OMEGA2D  — 

Maximum  depression  angle  in  beam  cone  (deg) 

DGRZD  — 

Angle  increment  (deg)  in  incident  fan  of  paths 
from  scattering  point  to  source  cylinder 

Record  6:  (Z(L),  C(L),  L 

=  1,  NSVP) 

Z  — 

Array  of  depths  (m)  for  sound  speed  profile 

C  — 

Array  of  corresponding  sound  speeds  (m/s) 
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